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Abstract 

The main goal of this paper is to assess the Hmits of vahdity, in the regime of low concentration 
and strong Coulomb coupling (high molecular charges), for a simple perturbative approximation 
to the radial distribution functions (RDF), based upon a low-density expansion of the potential of 
mean force and proposed to describe protein-protein interactions in a recent Small- Angle-Scattering 
(SAS) experimental study. A highly simplified Yukawa (screened Coulomb) model of monomers 
and dimers of a charged globular protein (/3-lactoglobulin) in solution is considered. We test the ac- 
curacy of the RDF approximation, as a necessary complementary part of the previous experimental 
investigation, by comparison with the fluid structure predicted by approximate integral equations 
and exact Monte Carlo (MC) simulations. In the MC calculations, an Ewald construction for 
Yukawa potentials has been used to take into account the long-range part of the interactions in the 
weakly screened cases. Our results confirm that the perturbative first-order approximation is valid 
for this system even at strong Coulomb coupling, provided that the screening is not too weak (i.e., 
for Debye length smaller than monomer radius). A comparison of the MC results with integral 
equation calculations shows that both the hypernetted-chain (HNC) and the Percus-Yevick (FY) 
closures have a satisfactory behavior under these regimes, with the HNC being superior throughout. 
The relevance of our findings for interpreting SAS results is also discussed. 

PACS numbers: 05.20.Jj, 61.10.Eq, 61.12.Ex 
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I. INTRODUCTION 



In spite of the large effort devoted in the last decades, a clear understanding of the 
interactions of macromolecules in solution is still far from being achieved . In particular, 
this is true in the case of globular proteins, which share with colloidal systems a number of 
common properties P|. 

^From the experimental point of view, there exist several biophysical techniques for ob- 
taining quantitative data on protein-protein interactions under physiologically relevant con- 
ditions. Small-angle scattering (SAS), for instance, is currently believed to provide very 
reliable information, under very different experimental conditions (pH, ionic strength, tem- 
perature, etc.). If the particle form factors are known, dividing the SAS intensity by the 
average form factor yields the experimental average structure factor, which is related to the 
radial distribution functions (RDF) gij{r) {i and j are species indexes). A recent experiment 
jj] reported Small-Angle- X-ray Scattering (SAXS) measurements on structural properties 
of a particular globular protein, the /3-lactoglobulin (/5LG), in acidic solutions (pH = 2.3), 
at several values of ionic strength in the range 7-507 mM. For this protein there is a clear 
evidence of a monomer-dimer equilibrium affected by the ionic strength of the solution 
P], and the authors of Ref. 0| were able to achieve a good fit of the experimental data 
by using a highly simplified "two-component macroion model" (mimicking monomers and 
dimers of /3LG), with effective forces represented by hard-sphere (HS) terms plus the repul- 
sive Yukawa (screened Coulomb) part of the well-known Derjaguin-Landau-Vervey-Overbeek 
(DLVO) potential Q]. One important novelty of that study, compared with previous ones, 
is the proposal of a relatively simple, improved approximation to the RDFs, suitable for 
best-fit programs and not restricted to the particular model but equally well applicable to 
different spherically-symmetric potentials. 

iFiom. the theoretical point of view, information on intermolecular forces can be extracted 
from the experimental average structure factor by comparison with a theoretical one, whose 
calculation requires the choice not only of an interaction model but also of a recipe for deriv- 
ing the RDFs from the intermolecular potentials. At present, the most accurate techniques 
for evaluating RDFs are the "exact" computer simulations - Monte Carlo (MC) and molec- 
ular dynamics (MD) - and the "approximate" integral equations (IE) from the statistical 
mechanical theory of classical fluids 0. Unfortunately, such complex methods can hardly 



3 



be included into a best-fit program for analyzing experimental data. In fact, MC or MD 
simulations require long computational times, and become difficult in regimes characteristic 
of globular proteins in solution (i.e., low concentration, high charges, asymmetry in size and 
charge among the components of the mixture). On the other hand, only for a very limited 
number of simple potentials and within an even more limited number of approximate "clo- 
sures", lEs of liquid theory admit analytical solutions, providing closed-form expressions 

cases, an iterative numerical procedure 
is necessary, and this poses a major drawback to any fitting scheme. Moreover, numerical 
solution of the IE closures tends to become unstable or does not converge in the region of 
our interest. 

In order to simplify the problem, most analyses of SAS data for highly dilute solutions 
employ the crude approximation of neglecting all intermolecular forces, assuming either 
large interparticle separations or weak interactions. In this case, gij{r) = 1, the average 
structure factor equals unity and the SAS intensity depends only upon the average form 
factor. A common first improvement over the previous choice then corresponds to ap- 
proximating the RDFs with their zero-density limit, given by the Boltzmann factor, i.e. 
gij{r) = exp [— /3</)jj(r)], where (pijir) is the pair potential, and (3 = (ksT)'^ the inverse 
of the thermal energy (absolute temperature multiplied by Boltzmann's constant). How- 
ever, this zero-density approximation becomes insufficient at moderate concentrations or in 
regimes of colloidal or protein solutions when electrostatic interactions are strong, i.e. at 
low ionic strength. 

Motivated by this scenario, Ref. [J] proposed a more accurate representation of gij{r) that 
takes into account, according to a perturbative scheme, terms up to the first-order in the 
density expansion of the potential of mean force, Wij{r) = — j3~^ In gij{r) P] (note that the 
iyjj(r)-expansion should not be confused with the RDF one, since these two expansions differ 
even at the first-order in density). While the satisfactory best-fit results of Ref. jj] seem 
to indicate that a first-order approximation to Wij{r) (Wl-approximation) is sufficiently 
accurate for low concentrations such as the experimental conditions under study, there is no 
way, a priori, to tell where this approximation breaks down, in the absence of some "exact 
results" to compare with. On the other hand, as the experimental conditions present in the 
analysis of Ref. 0| are fairly typical in the context of proteins in solution, we feel that it 
would be interesting to make such a comparison. Thus the main subject of the present paper. 
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which complements the methodological part of the work of Ref. ^, is not the proposal of 
a new potential model for jSLG, but a test of the Wl-approximation against more accurate 
(MC and IE) structural results. 

We perform MC simulations, at constant volume, temperature and total number of 
macroparticles, for the same HS-Yukawa-DLVO binary model, representing monomers and 
dimers of /9LG, investigated in Ref. P|. Various values of the screening parameter are con- 
sidered, and the MC results for gij{r) are compared with the corresponding ones predicted 
by the aforesaid Wl-approximation as well as by some commonly used lEs. In order to 
ensure always a good accuracy, the MC calculations are carried out with and without a 
suitable Ewald construction 11, which is expected to play a major role in the cases 
of strong long-range interactions (weak screening). Although the theoretical framework for 
the Ewald construction, well-known for unscreened Coulomb forces, has already been ex- 



12l |. this work represents, to the best of 



tended to Yukawa interactions in recent Refs. 
our knowledge, the first MC detailed analysis of its implementation and performance for the 
repulsive Yukawa case [l^. 

Our calculations allow a rather precise determination of the limits of validity for the Wl- 
expansion. They also show clearly the degree of reliability of some typical lEs under these 
frequently encountered, demanding, regimes. It is worthwhile stressing that our results are 
in fact rather general, as there exists a large variety of physical phenomena which can be 
described by Yukawa potentials jlj]. The existence of "exact" computer simulations for 
a binary model with these potentials would then prove to be useful within a much more 
general context than the one treated here. 



II. THE PROTEIN-PROTEIN INTERACTION POTENTIAL 

When mesoscopic (colloidal or protein) particles with ionizable surface groups are put 
into a microscopic polar solvent (like water), most of the charged surface groups dissociate 
into the solvent and form microscopic counterions, usually carrying one or two elementary 
charges. Consequently, the big particles acquire high charges of opposite sign and are called 
macroions or polyions. At equilibrium the counterions are located around the charged 
macroions, forming an electric double layer. The counterion distribution tends to screen 
the repulsions between macroions, which have charges of the same sign. The result is a 
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screened Coulomb (Yukawa) repulsion between macroions, which ensures the stability of 
the solution (charge-stabilization) with respect to a possible irreversible fiocculation. An 
important feature of such repulsions is that they can be tuned, by adding a suitable amount 
of a simple electrolyte to the solution. In fact, such a salt provides additional free microions 
(co-ions, with same charge sign as the macroions, as well as other counter-ions) j_which 
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increase the degree of screening and thus reduce the macroion-macroion repulsions 

A /5LG solution thus consists of many components: two different forms of macroions 
(protein monomers and dimers), counterions, colons and the solvent. At neutral pH, the 
structure of the (3LG protein is dimeric, while at acidic pH (a condition more similar to 
the physiological one) a partial dissociation into two monomers takes place. The monomer- 
dimer equilibrium, which determines the molar fractions of both macroion species, depends 
upon the ionic strength of the solution. At low ionic strength, the screening is weak and the 
electrostatic repulsions predominate over the attractive forces responsible for the formation 
of dimers; as a consequence, most of the macroions are monomers. On the contrary, at high 
ionic strength a strong screening reduces the monomer-monomer repulsions in such a way 
that a large fraction of dimers can form. 

As in Ref. 0|, we represent such a /?LG multicomponent solution at a highly simpli- 
fied, "primitive model", level of description, using an effective "two- component macroion 
model", which takes into account only protein particles Q, 3|- In fact, the solvent 
is regarded as a uniform dielectric continuum, all microions are treated as point-like parti- 
cles, and macroions (both monomers and dimers) are assumed to be charged hard spheres, 
with different diameters. The presence of both solvent and microions appears only in the 
macroion-macroion effective potentials. In the spirit of the DLVO theory p, we shall then 
describe the protein-protein interactions with the simple effective potential 

Mr) = <Pf{r) + <Pl,{r) (1) 

= 1, 2, with species 1 and 2 corresponding to monomers and dimers, respectively). Here 
the hard-sphere term accounts for excluded volume effects 

4 , (2) 

0, r > aij 

where aij = {cTi + o'j)/2 is the distance of closest approach between two macroparticles of 
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species i and j. On the other hand, the renormalized Yukawa term 



[r 



(3) 



"^'^^ ' e{l + KDai/2){l + KDaj/2) r 

represents an effective screened Coulomb repulsion between two isolated macroions in a sea 
of microions, and has the same Yukawa form as in the Debye-Hiickel theory of electrolytes, 
but with coupling coefficients of DLVO type j^. Here, e is the elementary charge, e the 
dielectric constant of the solvent, Zi the valency of species i, and kd is the inverse Debye 
screening length due only to microions, given by 



«:i. = \/^Y5^(/c + /.). (4) 

Na is the Avogadro number, and Jc = {l/2)ccZ^ denotes the ionic strength of the counterions 
originated from the ionization of the protein macromolecules (the molar concentration Cc of 
these counterions is related to the macroion concentrations through the electroneutrality 
condition, |Z,| = ci |Zi| + ca IZ2I), while /, = (1/2) c™ (^™)^ is the ionic strength 
of all microions (cations and anions) generated by added salts. Clearly, depends on 
temperature and represents an indication of the range of the screened Coulomb interactions, 
with kd ^ corresponding to pure Coulomb potentials, whereas k/j ^ 00 yields the 
opposite HS limit. While in real experiments is fixed by the chemical conditions of the 
solution (namely Ic and Ig), in this work we shall not consider Ig as an independent variable, 
but, in view of our methodological purpose, we shall regard k/jcti = C as an independent 
reduced screening parameter. 

A measure of the concentration of the two-macroion effective mixture will be given by 
the volume fraction 

TT ^ 

V = TT.Pi'^f (5) 

where pt is the partial number density of the z-th macroion species (point-like microions and 
solvent do not appear here). The definition of the model is then completed by providing one 
of the two molar fractions Xi = Pi/ p {i = 1,2), where p = J2iPi is the total density. 

Now, following partly Ref. jj], we add three remarks about some assumptions involved 
in the choice of the model potential. 

i) At first glance one might suspect that reducing dimers to equivalent spheres (with a 
volume twice as large as the monomer), i.e. neglecting the asymmetry of the dimer molecular 
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shape may seem a too drastic simplification. In order to clarify this point, it is to be stressed 
that in Refs. two different levels of description for the dimer were used in the two 

factors which contribute to the SAS intensity. The coherent scattering intensity /(g) was 
written as 

Iiq)^Y.M'^'^:i^)^Ms^M^ (6) 

where q is the magnitude of the scattering vector, -Fj(g) the angular average of the form factor 
of species i, and the Ashcroft-Langreth partial structure factors (for spherically-symmetric 
intermolecular potentials) are defined by 

S,,{q) = 5,, + 4n {p.p^f r r%^(rf-^^^dr (7) 

Jo qr 

in terms of the three-dimensional Fourier transform of hij{r) = gij{r) — 1. A very accurate 
procedure was used to calculate numerically both macroion form factors, Fi{q) and F2{q), 
from crystallographic data, taking into account, in particular, the exact elongated shape 
and structure of the dimer, i.e. its distribution of scattering matter 0, 0|. Thus the approx- 
imation of spherical dimers was restricted only to the calculation of Sij{q), which is related, 
through gij{r), to the intermolecular potentials. At low protein concentrations, the choice 
of spherically-symmetric hard-core potentials can indeed be justified. As in such regimes 
the average distance among particles is large, intermolecular forces are dominated by the 
long-range electrostatic interactions, whereas the details of the short-range repulsions (i.e. 
the excluded volume effects) are irrelevant. 

ii) Our potentials are purely repulsive. We have not included the attractive van der Waals 
)art of the DLVO potential for charged colloidal suspensions (the so-called Hamaker term 
^), as it has already been shown to be unnecessary for this system in previous work j^. 
The basic reason is that van der Waals attractions may be fully masked by the electrostatic 
repulsions when the latter are strong, and are also negligible for moderately charged particles 
with a diameter smaller then 50 nm Q]. Moreover, the Hamaker term diverges at contact, 
so that, to circumvent this singularity, the inclusion of the attractive term would require the 
addition of a Stern layer of counterions (with finite size) condensed on the macroion surface 

a. 

iii) Given that the specific protein forms dimers, it appears that the j3LG necessarily 
has a short-range monomer-monomer attraction (related to the surface groups), that causes 
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the aggregation into dimers and determines the monomer molar fraction xi. One expects 
this attractive term (possibly including hydrogen bonding) to be rather complex and non- 
spherically-symmetric. If such a contribution were clearly understood and easily tractable, 
one could start from a more fundamental viewpoint, choosing a model which considers only 
monomers and includes the aforementioned attraction into their pair potential. One could 
then monitor the dimerization fraction within this one- component system. However, this 
analysis may be a project on its own right and goes beyond the aims of the present study. 

More simply, in order to avoid poorly known and angular-dependent potentials, the au- 
thors of Refs. y, Ol adopted the viewpoint of using a binary (monomer-dimer) rather than a 
one-component model, and the required attraction was accounted only indirectly, by using 



a chemical association equilibrium to evaluate Xi 

While the dependence of xi upon the added salt (i.e. upon Ig) must be taken into account 
in any best-fit analysis with the binary model 0, S], in the present work, for the sake of 
simplicity, we shall consider xi as an independent parameter. Most of our calculations will 
be performed at equal molar fractions, xi = X2, but in the last part of the paper we shall 
also address the effect of changing the molar fractions. 

III. LOW-DENSITY EXPANSION OF THE MEAN FORCE POTENTIAL 

As discussed in the introduction, one of the most commonly used procedures to compute 
RDFs gij{r) for a given pair potential (f)ij{r) goes through the solution of the Ornstein- 
Zernike (OZ) lEs from the liquid state theory, within some approximate closure relation. 
This can typically be done only numerically, with the exception of few simple cases (for some 
potentials and peculiar closures), where the solution can be worked out analytically 

ote^at, for HS- Yukawa potentials, the OZ equations do admit analytical solution 
within the so-called "mean spherical approximation" (MSA), to be discussed 



Note thi 



further below. On the other hand, under the experimental regime which we are interested 



m 



namely low density and strong electrostatic repulsions (weak screening), the MSA 
is well known to display a serious drawback since RDFs may assume unphysical negative 
values close to contact distance (Xij, for particles i and j which repel each other. To overcome 
this shortcoming for repulsive Yukawa models, it would be possible to utilize an analytical 



"rescaled MSA" 



21 



22I ] (this possibility will not be investigated in the present paper) 
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or to resort to different closures. 

In general then, only numerical solutions are feasible, and thus IE algorithms can hardly 
be included into best-fit programs for the analysis of SAS results. 

The use of analytical solutions, or simple approximations requiring only a minor compu- 
tational effort, is clearly much more advantageous when fitting experimental data. This can 
be done by resorting to the following exact, albeit formal, relation 



g,, (r) = exp [-^3 (^)] , (8) 

- l3Wij (r) = -(3(t),j (r) + (r) , (9) 

where Wij (r) is the potential of mean force, which includes the direct pair potential (^ij (r) 
as well as —(3~^ujij{r), i.e. the indirect interaction between i and j due to their interactions 
with all remaining macroparticles of the fluid. In the density expansion of Wij (r) 

- m, {r) = Ins., (r) = (r) + 4f{r)p + 4f{r)p^ + ..., (10) 

the ..a.t powe. coefficients ^ 1, 2. . . .) can be computed by u.„g standard dia- 

grammatic techniques [9|| , which yield the results in terms of multi-dimensional integrals of 
products of Mayer functions 

^•(r)=exp[-/?0,,(r)]-l. (11) 

In the zero-density limit, Wij{r) vanishes and Qij (r) reduces to the Boltzmann factor, i.e. 

9ij {r) = exp [-/3(f)ij (r)] as p ^ 0, (12) 

which represents a 0*^-order (WO) approximation, frequently used in the analysis of ex- 
perimental scattering data. The WO- approximation avoids the problem of solving the OZ 
equations, but is largely inaccurate except, perhaps, at extremely low densities. We then 
consider the I'^^-order perturbative correction (Wl-approximation) jj] 

9ij (r) = exp [-(3(j)ij (r) + wj]^(r)p] . (13) 

By construction, this expression is never negative, thus overcoming the major drawback of 
MSA. The explicit expression of ujlj\r) reads 

^ff(r) = E^^7ff.(r) = / dr' (r') (|r - r'|) . (14) 

k k •' 
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The evaluation of the convolution integral 'yij\{r) is most easily carried out in bipolar coor- 
dinates. After an integration over angle variables 7iJl.(T') reduces to 

^Skir) = - dx [xU (x)] dy [yhM- (15) 

/ JO J\x—r\ 



Of course, the use of the Wl-approximation is not restricted to the model of this paper, 
and the proposed calculation scheme can be equallv well applied to different spherically- 
symmetric potentials. While it was shown in Ref. ^ how this first-order correction largely 
improves the fit of experimental scattering data, over the WO-one and under those exper- 
imental conditions, little could be said on the limits of validity of the Wl-approximation 
with respect to an (hypothetical) exact calculation. This is the reason why we tackle this 
task here by a comparison with MC simulations for a binary HS-Yukawa-DLVO system. 



IV. MC SIMULATIONS AND EWALD SUM FOR YUKAWA FLUIDS 

The difficulties involved in Monte Carlo calculations dealing with pure Coulomb potentials 
are well known [23]. It is now widely appreciated the usefulness of the so-called Ewald sum 
for long-range electrostatic interactions |lfl[ On the other hand, a similar construction 
for Yukawa potentials has appeared in the literature quite recently IjJ, • We now briefly 
recall the procedure detailed in Refs. |ll|,[3|. In order to keep notation as simple as possible, 
we shall restrict ourselves to simple Yukawa potentials, the extension to our actual potential 
(Eqs. [DEI) being obvious. The basic idea is to start with the total potential energy 

^ = 9 E ^"9/3— ' (16) 

where N is the total number of macroparticles, rap = Ifq, — r/3| and qa = Zae, qp = Zpe 
are the charges. This term is then split into a sum of two contributions, one evaluated 
in real space, while the other is calculated in momentum space on wave vectors given by 
k = 27rn/y^/^ {V is the volume of the system and n a unit vector of integer components). 
To this aim an auxiliary continuous Gaussian charge distribution 



\2\ 3/2 

^ ) e-''^' (17) 



TT / 
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is exploited. For A values such that the real space contribution is limited to particles in the 
basic simulation cell, the final result reads 



N 



u 



J2 '^o.qp 



erfc (Ar„^ + e'^^^'^i' + erfc (Ar„^ - e'''^'''^^ 



^0/3=1 
N 



^ 1 ^ 47r 



exp 



a/3: 



2A 



exp 



-Hp 

"Ia^ 



1^ 4A2 
erfc 



cos (k, 



a/3 ■ fa/? J 



(18) 



where in the first sum we exclude the terms with equal indexes and we have introduced the 
complementary error function 

2 



erfc(x) = —j= / dz e 



(19) 



The first two terms in Eq. (fTH|) represent the real and momentum space summations respec 



tively, while the last two contributions refer to the self-energy [Hi, 113 • In the limit 0, 



im spa c 

HQ. 



the above equation reduces to the Coulomb case [2^, as it should. Eq. (|TH|l contains A as 
an adjustable parameter, and we have performed a detailed analysis for its optimal choice, 
so that the original potential ()16|) is recovered for the range oi values of interest under 
the experimental conditions of Ref. ^, without using too many terms in the reciprocal 
space summation. Our results indicate A ~ 6.5/L (with L being the side length of the cubic 
simulation box) to be the optimal choice, which is of the same order of magnitude of the 
one typically used in the Coulomb case. 



V. INTEGRAL EQUATIONS 

Our next task is to test the performance of some lEs under the experimental conditions 
of Ref. j^. This will strengthen the usefulness of the Wl-approximation, in view of its 
simplicity compared to a typical IE calculation for a binary mixture. The OZ lEs of the 
iquid state theory for p-component mixtures with spherically-symmetric interactions read 



hij (r) = dj (r) 



pj^^i I dr'cKr') hi,{\v-v'\), (20) 
1=1 '' 
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and their solution can be accomplished only in the presence of an additional approximate 
relation (closure) between the direct correlation function (DCF) Cij (r) and the total corre- 
lation function hij (r) = — 1 (p = 2 in the present case). The most known among 
these approximations are 7[ 



1) The Percus-Yevick (PY) closure 



Cij [r) 



[1 + {r)] 



(21) 



where (r) = hij (r) - dj (r) . 

2) The Hypernetted Chain (HNC) closure 



(22) 



3) The Mean Spherical Approximation (MSA), much simpler than the above two, with the 
DCF being related only to the potential outside the core 



Cii {r) = -P(t>ij (r) r > 



I] -I 



(23) 



complemented by the condition of excluded volume, gij{r) = inside the hard cores 

Other possible more refined closures, which can be regarded as a combination of the 
above three, will be also briefly addressed in this work. 

VI. NUMERICAL RESULTS 



A /?LG-monomer is composed of 162 amino acid residues; 20 of these are basic, so that at 
pH = 2.3 the monomer is expected to be positively charged, with about 20 proton charges. 
In our calculations we fix all parameters close to their best-fit "experimental" values j^, 
ai = 40 A, a2 = ~ 50. 40 A , Zi = 20, Z2 = 40, T = 298.15 K and e = 78.5 (strictly 

speaking, in Ref. T = 293.15, cti = 38.30 A, and the ratio Z2/Z1 was about 1.8, since 
two of the 20 amino acids of the monomer are at the monomer-monomer interface in the 
dimer) . 



The packing fraction rj = 0.01 is a 



so very close to that determined from the experimental 



protein concentration (77 = 0.0096) J|. We then vary the dimensionless screening parameter 
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( = kdO'i in the range C ~ 1 — 10, roughly equivalent to the range of ionic strength Ig (from 
7 to 507 mM) examined in the aforesaid SAS measurements for jSLG [ where ( = 1.41 
when Is = 7 mM (weak screening, monomer molar fraction Xi = 0.85), and C = 9.08 when 
Is = 507 mM (strong screening, Xi = 0.05) ]. Note that an increase of ( has the effect of 
reducing not only the range of the HS-Yukawa-DLVO potentials but also their amplitudes, 
as described by Eq. Q. 

In order to obtain the Wl- approximation to the RDFs, we have evaluated all the convo- 
lution terms 7i|fc(r), given by Eq. (fT3j). at the grid points rj = iAr {i = 1, . . . , 500), with 
Ar = lA. At each value, the double integral Eq. ()15|l has been carried out numerically, 
by using the trapezoidal rule for both x— and y- integration. For the x-integration, we have 
chosen as upper limit the value x^a.^ = niax(xcut5 c"2 + f)^ with Xcnt = cr-i + ^'^I^d-, and as 
grid size Ax = Xcut/400. For the y-integration. Ay = Ax. 

The MC simulations have been performed at constant A^, V, T, with and without the 
Ewald procedure for a correct treatment of the long-range electrostatic interactions. Most 
calculations refer to a total number of particles A^ = 216, divided in monomers and dimers 
according to the fixed monomer molar fraction x\. Although the sample size may seem 
rather small with respect to present-day standards, one has to take into account that the 
Ewald construction takes a great computational effort with increasing A^. In any case, we 
have carried out some additional calculations with a larger number of particles in order 
to check for possible finite-size effects, and found no significant differences in the results. 
Hence, we shall use this value of A^ throughout, with one exception which will be described 
later on. The simulation starts from an appropriate lattice distribution of molecules. We 
have typically employed 10^ equilibration steps to eliminate any memory of the initial con- 
figuration artificially introduced into the fluid. Then 5 x 10^ additional steps have been used 
to collect sufficient information for the statistical averages required to calculate the RDFs. 

With the same parameters we have also solved the OZ i nteg ral equations numerically, by 
means of an efficient algorithm proposed by Labik et al. [2J| employing 1024 grid points, 
with a mesh size Ar = O.Olcri, and 20 basis functions. The PY, HNC and MSA closures 
have been employed. As expected, the MSA results (not shown in our Figures) poorly 
describe the MC data and exhibit the above-mentioned drawbacks of the MSA closure in 
regimes with strong coupling at high dilution jQ. We have explicitly checked that other, 
more sophisticated, approximations, such as the the Rogers- Young (RY) closure 2^ or the 
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Zerah-Hansen (HMSA) one [26[, which attempt to achieve thermodynamic consistency of 
compressibihty and virial pressures by interpolating between two of the above closures (PY- 
HNC and MSA-HNC, respectively), are of no use here, in such a thermodynamic consistency 
is never achieved, presumably because of the combined effect of low densities and strong 



long-range repulsions 



27|. 



Finally, both MC and IE calculations for gij{r) have been compared with the corre- 
sponding results from the first-order Wl-approximation, with the aim to assess the limits 
of validity where the expression given by Eq. (fTSj) can be safely exploited, under conditions 
typical of proteins in solution. As further elaborated below, we find that for values C ~ 2 
(i.e. ~ ^1/2) the Wl-approximation well describes the behavior of the RDFs. 

When C is large (in the range C ~ 5 — 10) the Yukawa interactions are strongly screened, 
and the RDFs essentially reduce to the typical HS ones, with the first maximum correspond- 
ing to the contact distance aij. 

Fig. Q] depicts the comparison between the MC results and the Wl-approximation for 
C = 3 (corresponding to a moderately weak screening) and xi = 0.5, that is when both 
monomers and dimers are present in equal measure. Note that these conditions are close to 
one of the experimental cases reported in Ref. Q], where Ig = 47 mM corresponds to ( = 2.8 
and Xi = 0.48. On the other hand, as the ionic strength Is is lowered from 507 mM to 7 
mM, the experimental system switches from a fluid almost completely made up of dimers 
(xi = 0.05) to one almost completely made up of monomers (xi = 0.85). This rather peculiar 
feature is specific of the f3LG and will also be considered further on. Here, however, our 
main aim is to test the Wl-approximation under the simple, symmetric, condition of equal 
molar fractions, since we already know, from Ref. 0|, that the first-order approximation 
well describes the (3LG experimental data, which display, in particular, a lowering in the 
scattering intensity at small angles, with a progressive development of an interference peak 
at low ionic strengths. In Fig. Q we also report the results from the HNC and PY lEs 
(solid and dotted lines), which are practically indistinguishable on the employed scale. It is 
apparent that in the case of Fig. [T]the Wl-RDFs gii{r), guir) and g22{r) are in excellent 
agreement with their MC, HNC and PY counterparts. Note that, for all three RDFs, gij{r) 
remains zero even in a region outside the hard-core, while the position of the peak lies at a 
distance larger than a^j, as a consequence of the strong Yukawa repulsions. 

A departure of the first-order Wl-approximation from the MC results can be observed 
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for smaller values of the screening parameter (, where higher-order terms in the density 
expansion of Wij (r), Eq. pO|) . begin to have a non-negligible effect. This is indicated in 
Fig. |21for the case C = 2, which corresponds to = cri/2, with the Debye screening length 
being equal to the monomer radius (among the experimental data of Ref. we find C = 2 
and Xi = 0.73 when Ig = 17 mM). Again the HNC and PY RDFs are nearly identical with 
each other and with MC data. On the other hand, the Wl-approximation predicts peaks 
nearly at the same positions as the PY and HNC closures, while its peak heights are slightly 
overestimated. However, the agreement between Wl and MC results can still be regarded 
as rather good. 

In regimes with weaker screening the discrepancies become more and more pronounced. 
The breakdown of all the considered approximations can be clearly appreciated in Fig. |21 
for C = 1 (note that the case with the weakest screening in Ref. ^ corresponds to = 7 
mM, C = 1.41 and Xi = 0.85). The Wl-results are not reported in this Figure, since they 
are way off from the MC data (with an overestimation of about a factor 2). On the other 
hand, even the results from the PY approximation are significantly displaced from the MC 
RDFs. The difference between the HNC and PY results is apparent, particularly for the 
latter, as expected. The PY approximation overestimates both the heights and positions 
of the peaks, compared to the HNC ones. Overall the PY approximation fails to describe 
the MC calculation for C < 2, whereas the HNC closure is consistently in good agreement 
with the MC data. Such a good performance of the HNC closure closely resembles the 
good agreement between HNC and MC, even at strong Coulomb coupling, for the one- 
component fluid of point charges (electron gas or plasma, with ( = 0) (OCP) in a uniform 
neutralizing background 



28[. However, the results for our binary model with screening at 
packing fraction rj = 0.01 can hardly be compared with the available MC simulations for one- 
component charged hard spheres (OCCS, with = 0) in a uniform neutralizing background, 
at ?7 = 0.3 -T- 0.4 [2^. Moreover, it is known the inadequacy of the HNC for high charges 
at low concentrations (for instance, in the dilute regime of 2-2 aqueous electrolytes 0, 0| , 
where bridge diagrams become non- negligible for like-charge RDFs). On the other hand, 
despite the large number of comparisons among PY, HNC and MC predictions carried out 
over the years, we are not aware of a similar detailed RDF investigation under regimes 
characteristic of globular proteins in solution, for HS-Yukawa-DLVO binary models. 

Next we consider the effect of taking into proper account the long-range nature of the 
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interactions (in the weakly screened case) with the use of the Ewald construction. This 
is illustrated in Figs. |3 and where the RDFs computed with and without the Ewald 
construction are compared at ( = 1 and ( = 0.25, respectively. Clearly, very little difference 
is detected between these two calculations when C = 1 (and when ( = 0.5, not shown). 
We find that the the presence of the Ewald construction begins to be important for very 
low values of the screening parameter (C < 0.25, i.e. k]j^ > 4cri), as shown in Fig. |31 
Supplementary calculations, not reported here, confirm that this is true even for lower 
values of protein charges, that is for weaker Coulomb coupling. 

Finally, we consider the effect of varying the molar fractions. While the exact condi- 
tions reported in the (3LG experiment pose a very difficult challenge to an accurate MC 
calculation in view of the particular combination of strong asymmetry and repulsions, we 
can nevertheless easily account for the general trend. This is depicted in Fig. where we 
have assumed ( = 2 and Xi = 0.75, in closer analogy with a (3LG experimental case, ( = 2 
and Xi = 0.73 when = 17 mM. It is apparent how the performance of the first-order 
Wl-approximation is comparable to the corresponding symmetric case, C = 2 and xi = 0.5. 

Fig. [prefers to the asymmetric case with the weakest screening in Ref. [J], i.e. ( = lAl 
and Xi = 0.85 (corresponding to the lowest value of ionic strength, = 7 mM). Again, the 
HNC and PY results are in good agreement with the MC ones, and even the performance of 
the Wl-approximation can be regarded as acceptable, in agreement with the results of Ref. 
0]. We note that, in view of the low molar fraction of species 2 (dimers), the results of Fig. 
[prefer to a higher number of particles (A^ = 512). 



VII. CONCLUSIVE REMARKS 



This work represents a necessary verification of the best-fit analysis of SAS experimental 
data, for solutions of /3-lactoglobulin, presented in Ref. [J]. In the present paper we have 
assessed the limits of validity of the Wl-approximation, exploited in that work to calculate, 
in a simple way, the RDFs in regimes typical of a large class of globular proteins in solution, 
that is low concentrations and high macroion charges. This task has been accomplished by 
considering the same highly simplified model proposed in Ref. [J] (i.e. a binary mixture of 
monomers and dimers of the protein, with HS-Yukawa-DLVO effective potentials), and com- 
paring the corresponding gij{r) obtained by three different methods: the first-order density 
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expansion of the potential of mean force (Wl-approximation), "exact" MC simulations and 
approximate lEs. All results reported here refer to 77 = 0.01 and high macroion charges, 
Zx = 20 and Z2 = 40. For the MC simulations we have implemented an Ewald construction 
for Yukawa potentials, which ensures a proper treatment of the long-range part of the inter- 
actions, and we have tested its relevance as a function of the screening parameter (. In the 
IE calculations simple closures (PY, HNC, and MSA) as well as more elaborated ones (RY 
and HMSA) have been considered. 

We can summarize the obtained results as follows. 

i) The first-order Wl-approximation can be considered reliable in regimes with low con- 

centration [rj = 0.01) even for strong Coulomb coupling (up to charges of 10 ^ 20e on 
macroions with diameters of 40 ^50 A) , provided that the screening is strong enough, 
i.e. when C ^ 2 or, equivalently, ~ cri/2 (Debye length smaller than monomer 
radius). This finding demonstrates that the previous usage of the Wl-approximation 
in Ref. jj] was fully legitimate, for all considered cases including those with the lowest 
ionic strength (Jg = 7 mM, Xi = 0.85, ( = 1.41), which lies near the borderline of 
the reliability region. For weaker screening (lower values of C or larger k^^) at least 
second-order terms in the density expansion should be taken into account. However, 
the resulting W2- approximation would require a much higher computational effort 
and thus could not be conveniently included into a best-fit program for analyzing SAS 
experimental data. 

ii) In the MC simulations the Ewald construction for Yukawa potentials starts to be im- 

portant for weak screening corresponding to C ~ 0.25 (k^^ ^ 4cri), and this is true 
even for lower values of the protein charges. 

iii) Both the HNC and PY lEs yield sufficiently accurate values of the RDFs, as long as 

C ^ 2. For lower values of ( HNC is still accurate, whereas PY starts to deviate as 
expected. The MSA predictions, on the other hand, are very poor even in those regimes 
where the Wl-approximation can be considered reliable. Under these conditions both 
the RY and HMSA closures are found not to achieve thermodynamic consistency 
between compressibility and virial pressures. 

iv) The sufficient accuracy of the Wl-approximation in the regimes of our interest (tested 
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in this paper against "exact" MC results), together with its success (shown in Ref. 
j^) in reproducing the main features of the experimental SAS intensity curves for the 
examined /5LG solutions, confirm the good performance of the highly idealized two- 
macroion model, which includes spherically-symmetric HS-Yukawa-DLVO repulsions, 
a monomer-dimer chemical equilibrium, and the "exact" form factors, evaluated by 
taking into account the real non-spherical structure of the dimer. 

Clearly, all complex characteristics of the interactions between globular proteins cannot 



Q 



and here. We have 



be explained by the "primitive" level of description adopted in Ref. 
followed the generally accepted philosophy of exploiting the simplest possible description of 
the system, which yet can provide useful information on the basic underlying interaction 
mechanism. The determination of the "true" protein-protein potentials thus remains an 
open problem. 

Our choice of purely repulsive interactions illustrates the minimal assumptions allowing 
a satisfactory reproduction of the SAS data for /3LG. In many studies on colloidal or pro- 
tein solutions, satisfactory results were obtained from very simplified models. The use of 
sophisticated potentials, with a large number of different contributions, is often unnecessary 
at the first stages. Moreover, a high level of description for potentials would be in striking 
contrast with the poor level of approximation to the RDFs (WO-approximation) commonly 
adopted in many analyses of experimental data. 

As regards the approximation of spherical symmetry, used for the protein-protein interac- 
tions (but not in the calculation of the form factors), we remark that it represent a common 
simplifying choice. In particular, it is worth recalling a very recent study by Pellicane et 
al. jSj], which reports evidence that the phase diagram of prototype globular protein so- 
lutions (lysozyme and 7-crystallin in water and added salt) can be reasonably reproduced 
by a spherically-symmetric representation of macromolecular interactions. These authors 
employed a HS-Yukawa-DLVO one-component potential, including the Hamaker attractive 
part. 

Evidently, in addition to the molecular granularity of the solvent and the finite sizes of 
all microions, a highly refined model description of protein solutions should embody the 
asymmetry of the molecular shape as well as the heterogeneity of the macroion surface 
charge distribution. The presence of different charged surface groups may produce "charge 
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patches" that have a sign opposite to that of the net macroion charge. The importance 
of non-spherically-symmetric models with an inhomogeneous distribution of positively and 
negatively charged groups was recently investigated in a MC study on the electrostatic 
complexation of flexible polyelectrolytes with a-lactalbumin and /5-lactoglobulin 33^ . 

As a final remark to the present paper, it is worth pointing out that we are not aware 
of any previous investigations of this type within the HS-Yukawa-DLVO binary model and 
in regimes typical of globular proteins in solution. Our results and the methodological ap- 
proach based upon the Wl- approximation are expected to be useful in the analysis of SAS 
experiments. It would be rather interesting to pursue a similar study on the thermody- 
namic predictions of the first-order approximation. This could be easily carried out, as all 
thermodynamic quantities can be inferred either directly or through the knowledge of the 
RDFs. Another interesting issue, within the present framework, involves an increase of the 
asymmetry between the two considered molecular sizes, which is known to lead to possible 
depletion effects js^. We plan to perform such investigations in a future publication. 
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FIG. 1: Partial correlation functions gii{r), gi2{f) and g22if') (in order from bottom to top) as a 
function of the rescaled distance r/ai for C = 3 and xi = 0.5. Circles correspond to MC calcula- 
tions, full lines to HNC, dotted lines to PY, and dashed lines to the first-order Wl-approximation. 
Here and in the following the components 12 and 22 have been shifted upwards by one and two 
units, respectively. 
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FIG.2 Giacometti et al 




FIG. 2: Same as above with C = 2 and xi = 0.5. 
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FIG.3 Giacometti et al 




FIG. 3: Comparison of the results from HNC, PY and MC in the calculation of the partial radial 
distributions functions for C = 1 ^^iid = 0.5. The first-order Wl-approximation is not depicted 
as it overshoots the MC results roughly by a factor 2. 
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FIG.4 Giacometti et al 




FIG. 4: Partial correlation functions gii{r), gi2{r) and g22('") (in order from bottom to top) as 
a function of the rescaled distance r/ai, as computed with (circles) and without (solid line) the 
Edwald construction, for C = 1 = 
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FIG.5 Giacometti et al 




FIG. 5: Same as above with C = 0.25 and xi = 0.5. Note that the scale has been changed with 
respect to previous figures. Accordingly, here components 12 and 22 have been shifted upward by 
2 and 4 units, respectively. 



27 



FIG.6 Giacometti et al 




FIG.7 Giacometti et al 
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FIG. 7: Asymmetric case with (" = 1.41, xi = 0.85, corresponding to the lowest value of ionic 



strength Is = 7 mM (i.e., the weakest screening) investigated in Ref. 
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